.
.
Referências:
ggplotly(
ggplot(renner, aes(x = ref.date, y = ret.adjusted.prices)) +
geom_line() +
labs(x = 'data', y = 'retorno (em %)')
)
# guardando série de retorno
rt <- renner$ret.adjusted.prices
acf(rt)
pacf(rt)
Os gráficos indicam que a śerie \(r_t\) não possui autocorrelação.
Porém, pelo gráfico da série, podemos ver que a variância não parece ser constante ao longo do tempo, pois em alguns dias seguidos o retorno é muito alto ou muito baixo.
As séries correlacionadas podem ser estudadas usando o valor absoluto ou o valor quadrático. Se os retornos diários são i.i.d., então os retornos diários em valor absoluto (ou ao quadrado) também serão.
par(mfrow = c(2, 1))
acf(abs(rt), main = 'retorno em valor absoluto')
acf(rt ^ 2, main = 'retorno quadrático')
par(mfrow = c(1, 1))
par(mfrow = c(2, 1))
pacf(abs(rt), main = 'retorno em valor absoluto')
pacf(rt ^ 2, main = 'retorno quadrático')
par(mfrow = c(1, 1))
Como podemos ver, as séries em valor absoluto e ao quadrado possuem autocorrelação, contradizendo a hipótese (para a série original \(r_t\)) que não existia autocorrelação.
Os modelos ARCH, ou modelos auto-regressivos com heterocedasticidade condicional, foram introduzidos por Engle (1982), com o objetivo de estimar a variância da inflação. A expectativa aqui é que o retorno \(X_t\) seja não-correlacionado serialmente, mas a volatilidade (variância condicional) depende de retornos passados por meio de uma função quadrática.
Um modelo \(ARCH(p)\) é definido por
\(X_t = \sqrt{h_t}\epsilon_t\),
\(h_t = \alpha_0 + \alpha_1X_{t-1}^2+\ldots+ \alpha_rX_{t-p}^2\),
onde \(\epsilon_t\) é uma sequência de variáveis aleatórias independentes e identicamente distribuídas (i.i.d.) com média zero e variância um, \(\alpha_0 > 0\), \(\alpha_1 \geq 0\), \(i > 0\).
Na prática, usualmente supomos \(\epsilon_t\) ~ \(N(0,1)\) ou \(\epsilon_t\) ~ \(t_v\) (distribuição \(t\) de Student com \(v\) graus de liberdade).
\[ E(X_t) = 0 \]
\[ Var(X_t) = \frac{\alpha_0}{1-\alpha_1} \]
\[ Cov(X_{t+k}X_t|F_{t+k-1})=0 \]
O modelo ARCH na verdade é um modelo de regressão, em que a variável resposta é a volatilidade condicional e os lags passados de \(r_t ^ 2\) são as covariáveis.
O modelo \(ARCH(1)\) considera somente o retorno quadrático no instante \(t\) como função do retorno quadrático anterior, ou seja, \(r_{t-1} ^ 2\). Portanto num modelo \(ARCH(q)\), considera-se os \(t - q\) retornos quadráticos, e não somente o anterior.
Uma generalização dos modelos \(ARCH\) foi sugerida por Bollerslev (1986), chamado modelo \(GARCH\) (generalized \(ARCH\)).
Um modelo \(GARCH(p,q)\) é definido por
\(X_t = \sqrt{h_t}\epsilon_t\),
\(h_t = \alpha_0 + \sum_{\substack {i=1}}^{r} \alpha_iX_{\substack {t-i}}^{2} + \sum_{\substack {j=1}}^{s} \beta_jX_{\substack {t-j}}\) ,
em que \(\epsilon_t\) i.i.d., \(\alpha_0 > 0\), \(\alpha_i \geq 0\), \(\beta_j \geq 0\), \(\sum_{\substack {i=1}}^{q}(\alpha_i + \beta_i) <1\), \(r = max(p,q)\).
No modelo \(GARCH(p, q)\), além de considerarmos os retornos quadráticos passados, também leva-se em consideração as variâncias condicionais passadas.
O modelo \(GARCH(p, q)\) para a série original de retornos \(r_t\) implica que a série quadrática \(r_t ^ 2\) é modelada por um \(ARMA(max\{p, q\}, q)\). Então podemos usar todas técnicas de identificação que aprendemos até agora para modelar \(r_t^2\) como um modelo \(ARMA(p, q)\):
# deixando de fora 10 valores
rt_treino <- rt[1:(length(rt) - 10)]
# ajustando modelo ARMA
mod <- auto.arima(rt_treino ^ 2, d = 0)
summary(mod)
## Series: rt_treino^2
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.9838 -0.9677 3.9145
## s.e. 0.0127 0.0174 0.3853
##
## sigma^2 estimated as 47.44: log likelihood=-4127.46
## AIC=8262.93 AICc=8262.96 BIC=8283.4
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.009930523 6.879247 4.040285 -Inf Inf 0.7831633 0.01974922
O método automático auto.arima sugere um \(ARMA(1, 1)\) para ajustar o retorno quadrático \(r_t ^ 2\).
checkresiduals(mod, plot = F)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 9.3029, df = 7, p-value = 0.2316
##
## Model df: 3. Total lags used: 10
checkresiduals(mod, plot = F, lag = 12)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 9.8054, df = 9, p-value = 0.3665
##
## Model df: 3. Total lags used: 12
checkresiduals(mod, plot = F, lag = 16)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 17.202, df = 13, p-value = 0.1902
##
## Model df: 3. Total lags used: 16
checkresiduals(mod, plot = F, lag = 20)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 21.009, df = 17, p-value = 0.2259
##
## Model df: 3. Total lags used: 20
checkresiduals(mod, plot = F, lag = 24)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 23.412, df = 21, p-value = 0.3224
##
## Model df: 3. Total lags used: 24
Para diferentes lags, o teste de Ljung-Box não rejeita a hipótese nula de que os resíduos do modelo \(ARMA(1, 1)\) são ruído branco. Através desta identificação, podemos ajustar um \(GARCH(1, 1)\).
g_model <- garch(rt_treino, order = c(1, 1), trace = F)
summary(g_model)
##
## Call:
## garch(x = rt_treino, order = c(1, 1), trace = F)
##
## Model:
## GARCH(1,1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.14480 -0.55946 0.02364 0.69758 5.45180
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 0.08554 0.07141 1.198 0.2310
## a1 0.01682 0.00765 2.199 0.0279 *
## b1 0.96151 0.02432 39.541 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 78.439, df = 2, p-value < 2.2e-16
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.047388, df = 1, p-value = 0.8277
O livro Time Series Analysis (with Applications in R) - Jonathan D. Cryer, Kung-Sik Chan sugere que se o modelo \(GARCH(p, q)\), foi escolhido adequadamente, então o valor absoluto, ou quadrado dos resíduos, padronizados do modelo não devem possuir pontos significativos para o gráfico de ACF:
par(mfrow = c(2, 1))
acf(abs(residuals(g_model)), na.action = na.pass)
acf(residuals(g_model) ^ 2, na.action = na.pass)
par(mfrow = c(1, 1))
Os gráficos sugerem que o modelo \(GARCH(1, 1)\) é adequado aos dados de retorno de ações da Renner.
\(X_t = \sqrt{h_t}\epsilon_t\),
\(h_t = 0.08554 + \sum_{\substack {i=1}}^{q} 0.01682X_{\substack {t-i}}^{2} + \sum_{\substack {j=1}}^{s} 0.96151X_{\substack {t-j}}\) ,
em que \(\epsilon_t\) i.i.d., \(q = max(p,q)\).
# eleva ao quadrado porque o modelo ajusta o desvio padrão condicional, e não a variância condicional
ggplotly(
tibble(t = renner$ref.date[1:1233],
observado = rt_treino,
GARCH_pos = fitted(g_model)[, 1] ^ 2,
GARCH_neg = -fitted(g_model)[, 2] ^ 2) %>%
ggplot(aes(x = t)) +
geom_line(aes(y = observado), color = 'grey50') +
geom_line(aes(y = GARCH_pos), color = 'red') +
geom_line(aes(y = GARCH_neg), color = 'red4')
)
A predição retorna o desvio padrão condicional, por isso eleva-se ao quadrado.
# dados de teste
teste <- rt[1234:1243]
pred <- predict(g_model, newdata = teste)[, 1] ^ 2
tab_pred <- tibble(
t = renner$ref.date[1:1243],
observado = rt[1:1243],
predito_pos = c(rep(NA, 1233), pred),
predito_neg = c(rep(NA, 1233), -pred)
)
ggplotly(
tab_pred[1100:1243, ] %>%
ggplot(aes(x = t)) +
geom_line(aes(y = observado), color = 'grey50') +
geom_line(aes(y = predito_pos), color = 'red') +
geom_line(aes(y = predito_neg), color = 'red4') +
geom_hline(yintercept = 0, linetype = 'dashed')
)
rmse <- function(y, y_hat) {
sqrt(sum((y - y_hat) ^ 2) / length(y))
}
cat('RMSE para dados de treino:\n')
## RMSE para dados de treino:
rmse(rt[1235:1243], predict(g_model, newdata = rt[1234:1243])[, 1][2:10])
## [1] 1.773664
cat('RMSE para dados de teste:\n')
## RMSE para dados de teste:
rmse(teste[2:10], pred[2:10])
## [1] 3.424494